'''
WARNING:
  Preserve Python3 compatibility
  print()
'''


'''
Jorg Stelling model from
A synthetic tunable mammalian oscillator
'''

import sys
import math
from numpy import *
from scipy import *
from bricks import *

BA_BISTABLE     = 0
BA_MONOSTABLEMCT = 1
BA_MONOSTABLENEP = 2
BA_AMBIGOUS     = 3
BA_NOT_VALID    = -1

class BistableFeedbackSwitch:
  '''
  ARMUMENTS
    TIME : simulation time (minutes)
  '''
  gDEBUG = True
  gTIME = 10000
  gG1 = 200.0
  gG2 = 200.0
  gG3 = 200.0
  gG4 = 200.0
  gG5 = 200.0
  gG6 = 200.0
  gpc_start = [];
  gpc_half_life = 0;
  gpc_duration = []
  grm_start = [];
  grm_half_life = 0;
  grm_duration = []
  greq = 0
  gthresh = 1e2
  ghk = []
  

  '''
  MODEL PARAMETERS

  G1: [57][59][pmin][TAL97:KRAB][NEPTUN][TAL57:VP16]
  G2: [59][97][pmin][TAL57:KRAB][mCITRIN][TAL97:VP16]
  G3: [PIP][p(PIR)ON][TAL97:KRAB][TAL57:VP16]
  G4: [12xZF-BS][p-IL2][TAL57:KRAB][TAL97:VP16]
  G5: [pCMV][PIP-KRAB][HetAct]
  '''
  gZERO = 1e-4

  


  def __init__(self, dbg, t, g1, g2, g3,g4,g5, req, pcs, pchl, pcd, rms, rmhl, rmd, hk, thr):
    self.gDEBUG = dbg
    self.gTIME = t
    self.gG1 = g1
    self.gG2 = g2
    self.gG3 = g3
    self.gG4 = g4
    self.gG5 = g5
    self.gpc_start = pcs;
    self.gpc_half_life = pchl;
    self.grm_start = rms;
    self.grm_half_life = rmhl;
    self.gthresh = thr;
    self.gpc_duration = pcd;
    self.grm_duration = rmd;
    self.greq = req;
    self.ghk = hk;




  def plotResults(self,TIME, T97KRABm, T97KRABp, T57KRABm, T57KRABp, T97VP16m, T97VP16p, T57VP16m, T57VP16p,  PBSm, PBSp, PBSa, HACTm, HACTp, HACTa, PC, RM, NEPm, NEPp, MCTm, MCTp, BIFNEPp, BIFNEPt, BIFMCTp, BIFMCTt):
    if sys.version_info[0] > 2:
      return

    from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
    from matplotlib.figure import Figure
    from matplotlib.dates import DateFormatter
    from random import random
    
    NONCESIZE = 1000;
    NONCE = int(random()*NONCESIZE)

    if self.gDEBUG:
      f0 = figure()
      f1 = figure()
      f2 = figure()
      f3 = figure()
      f4 = figure()
      f5 = figure()
      
    else:
      f0 = Figure()
      f1 = Figure()
      f2 = Figure()
      f3 = Figure()
      f4 = Figure()
      f5 = Figure()
      
    f0.set_facecolor('white')
    s0=f0.add_subplot(111)
    s0.plot(TIME, HACTm, linestyle='dashed', label='HACT mRNA', color='blue')
    s0.plot(TIME, HACTp, label='HACT', color='blue')
    s0.plot(TIME, PBSm, linestyle='dashed', label='PIP mRNA', color='green')
    s0.plot(TIME, PBSp, label='PIP', color='green')
    s0.set_xlabel('Time [min]')
    s0.set_ylabel('Concentration [nM]')
    s0.legend(loc="upper right")
    lg = s0.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)

    f1.set_facecolor('white')
    s1=f1.add_subplot(111)
    s1.plot(TIME, PBSa, label='PIP-active', color='blue')
    s1.plot(TIME, HACTa, label='HACT-active', color='red')
    s1.plot(TIME, PC, label='Pristinamycin', color='blue', linestyle='dotted')  
    s1.plot(TIME, RM, label='Rapamycin', color='red', linestyle='dotted')
    s1.set_xlabel('Time [min]')
    s1.set_ylabel('Concentration [nM]')
    s1.legend(loc="upper right")
    lg = s1.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
     
    f2.set_facecolor('white')
    s2=f2.add_subplot(111)
    s2.plot(TIME, T97KRABm, linestyle='dashed', label='TAL97:KRAB mRNA', color='blue')
    s2.plot(TIME, T97KRABp, label='TAL97:KRAB', color='blue')
    s2.plot(TIME, T57KRABm, linestyle='dashed', label='TAL57:KRAB mRNA', color='green')
    s2.plot(TIME, T57KRABp, label='TAL57:KRAB', color='green')
    s2.set_xlabel('Time [min]')
    s2.set_ylabel('Concentration [nM]')
    s2.legend(loc="upper right")
    lg = s2.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    
    f3.set_facecolor('white')
    s3=f3.add_subplot(111)
    s3.plot(TIME, T97VP16m, linestyle='dashed', label='TAL97:VP16 mRNA', color='green')
    s3.plot(TIME, T97VP16p, label='TAL97:VP16', color='green')
    s3.plot(TIME, T57VP16m, linestyle='dashed', label='TAL57:VP16 mRNA', color='blue')
    s3.plot(TIME, T57VP16p, label='TAL57:VP16', color='blue')
    s3.set_xlabel('Time [min]')
    s3.set_ylabel('Concentration [nM]')
    s3.legend(loc="upper right")
    lg = s3.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    f4.set_facecolor('white')    
    s4=f4.add_subplot(111)
    s4.plot(TIME, NEPm, linestyle='dashed', label='NEPTUN mRNA', color='blue')
    s4.plot(TIME, NEPp, label='NEPTUN', color='blue')
    s4.plot(TIME, MCTm, linestyle='dashed', label='mCITRIN mRNA', color='green')
    s4.plot(TIME, MCTp, label='mCITRIN', color='green')
    s4.set_xlabel('Time [min]')
    s4.set_ylabel('Concentration [nM]')
    lg = s4.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)

    
    s5=f5.add_subplot(111)
    s5.plot(MCTp, NEPp, 'c.', label="Phase")
    s5.plot(BIFMCTp,BIFNEPp, 'rx', label="Stable state")
    s5.plot([self.gthresh for i in range(10000)],linspace(self.gZERO,1e4,10000), label="Threshold State #1", linestyle='dashed', color='gray');
    s5.plot(linspace(self.gZERO,1e4,10000), [self.gthresh for i in range(10000)], label="Threshold State #2", linestyle='dashed', color='gray');
    s5.loglog()
    s5.set_xlabel('mCITRIN')
    s5.set_ylabel('NEPTUN')
    s5.legend(loc="upper right")
    f5.set_facecolor('white')   
    lg = s5.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    
    if self.gDEBUG:
      path = './out/'+self.__class__.__name__
      paths = [path+'_s0_'+str(NONCE)+'.png', path+'_s1_'+str(NONCE)+'.png', path+'_s2_'+str(NONCE)+'.png', path+'_s3_'+str(NONCE)+'.png', path+'_s4_'+str(NONCE)+'.png', path+'_s5_'+str(NONCE)+'.png']
      f0.savefig(paths[0], format ='png')
      f1.savefig(paths[1], format ='png')
      f2.savefig(paths[2], format ='png')
      f3.savefig(paths[3], format ='png')
      f4.savefig(paths[4], format ='png')
      f5.savefig(paths[5], format ='png')
    
    else:
      path = '/var/www/imodel/img/'+self.__class__.__name__
      paths = [path+'_s0_'+str(NONCE)+'.png', path+'_s1_'+str(NONCE)+'.png', path+'_s2_'+str(NONCE)+'.png', path+'_s3_'+str(NONCE)+'.png', path+'_s4_'+str(NONCE)+'.png', path+'_s5_'+str(NONCE)+'.png']
      canvas=FigureCanvas(f0)
      f0.savefig(paths[0], format ='png')
      canvas=FigureCanvas(f1)
      f1.savefig(paths[1], format ='png')
      canvas=FigureCanvas(f2)
      f2.savefig(paths[2], format ='png')
      canvas=FigureCanvas(f3)
      f3.savefig(paths[3], format ='png')
      canvas=FigureCanvas(f4)
      f4.savefig(paths[4], format ='png')
      canvas=FigureCanvas(f5)
      f5.savefig(paths[5], format ='png')
    
    return paths


  def bifurcationAnalysis(self,BIFNEPp, BIFMCTp, NEPp, MCTp):
    
    
    BA_BISTABLE     = 0
    BA_MONOSTABLEMCT = 1
    BA_MONOSTABLENEP = 2
    BA_AMBIGOUS     = 3
    BA_NOT_VALID    = -1
    
    
    stateMCT = False
    stateNEP = False
    ret     = BA_NOT_VALID

    #Discover ambigous states
    for i in range(len(BIFMCTp)):
      a = (BIFMCTp[i] >= self.gthresh and BIFNEPp[i] < self.gthresh)
      b = (BIFMCTp[i] < self.gthresh and BIFNEPp[i] >= self.gthresh)
      if not (a or b):
      #if a == b:
        ret = BA_AMBIGOUS
        break;
      if BIFMCTp[i] >= self.gthresh:
        stateMCT = True
      if BIFNEPp[i] >= self.gthresh:
        stateNEP = True  

    if not ret == BA_AMBIGOUS:
      if stateMCT and stateNEP:
        ret = BA_BISTABLE
      elif stateMCT:
        ret = BA_MONOSTABLEMCT
      elif stateNEP:
        ret = BA_MONOSTABLENEP
      
    return ret

  
  #################
  # SIMULATION PART
  #################
  def run(self):
    TIME = range(self.gTIME)
    G1 = self.gG1
    G2 = self.gG2
    G3 = self.gG3
    G4 = self.gG4
    G5 = self.gG5


    #Experimental TAL Effector rates
    #RLA Measurements
    if self.gDEBUG:
      #tal97krab = Effector('TAL97:KRAB', 'dat/tal97krab.dat')
      #tal57krab = Effector('TAL57:KRAB', 'dat/tal57krab.dat')
      #tal97vp16 = Effector('TAL97:KRAB', 'dat/tal95vp16_4.dat')
      #tal57vp16 = Effector('TAL57:KRAB', 'dat/tal95vp16_4.dat')
      tal97krab = Effector('TAL97:KRAB', 'dat/talAkrab.dat')
      tal57krab = Effector('TAL57:KRAB', 'dat/talBkrab.dat')
      tal97vp16 = Effector('TAL97:KRAB', 'dat/talAvp16.dat')
      tal57vp16 = Effector('TAL57:KRAB', 'dat/talBvp16.dat')
      
    else:
      #tal97krab = Effector('TAL97:KRAB', '/home/root/cellulator-mirror/dat/tal97krab.dat')
      #tal57krab = Effector('TAL57:KRAB', '/home/root/cellulator-mirror/dat/tal57krab.dat')
      #tal97vp16 = Effector('TAL97:KRAB', '/home/root/cellulator-mirror/dat/tal95vp16_4.dat')
      #tal57vp16 = Effector('TAL57:KRAB', '/home/root/cellulator-mirror/dat/tal95vp16_4.dat')
      tal97krab = Effector('TAL97:KRAB', '/home/root/cellulator-mirror/dat/talAkrab.dat')
      tal57krab = Effector('TAL57:KRAB', '/home/root/cellulator-mirror/dat/talBkrab.dat')
      tal97vp16 = Effector('TAL97:KRAB', '/home/root/cellulator-mirror/dat/talAvp16.dat')
      tal57vp16 = Effector('TAL57:KRAB', '/home/root/cellulator-mirror/dat/talBvp16.dat')
      
    tal97krab.plotTransferFunction()
    tal57krab.plotTransferFunction()
    show()
    quit()
    
    #RATIOS
    r97krab = (G3+G1)/(G2)
    r57krab = (G2+G4)/(G1)
    r97vp16 = (G2+G4)/(G2)
    r57vp16 = (G3+G1)/(G1)

    #TAL:KRAB REPRESSORS
    t57max = self.gthresh;
    t97max = self.gthresh;
    hk97 = self.ghk[2]
    hk57 = self.ghk[3]
    if self.greq:
      #q = (tal97krab.getTransferRate2(r97krab)+tal57krab.getTransferRate2(r57krab))/2.0
      q97 = 0.1
      q57 = 0.1
    else:
      q97 = tal97krab.getTransferRate2(r97krab)
      q57 = tal57krab.getTransferRate2(r57krab)
    #print q97
    #print q57
    k57 = 1.0/(1/(q57)-1);
    k97 = 1.0/(1/(q97)-1);
    
    #TAL:VP16 ACIVATORS
    AMAX = 1383
    a97 = tal97vp16.getTransferRate2(r97vp16)
    a57 = tal57vp16.getTransferRate2(r57vp16)
    a97 = int(a97 > AMAX) * AMAX
    a57 = int(a57 > AMAX) * AMAX
    ha57 = hk57
    ha97 = hk97
    aMIN = AMAX
    Ka97 = self.gthresh
    Ka57 = self.gthresh
    
    '''
    tal97krab.plotTransferFunction()
    tal57krab.plotTransferFunction()
    tal97vp16.plotTransferFunction()
    tal57vp16.plotTransferFunction()
    print q97
    print q57
    print a97
    print a57
    show()
    quit()
    '''

    #Inducer affinity
    Kpc = 1.0
    Krm = 1.0

    #Maximal transcription rates (mRNA/min)
    #apply pCMV, pMin
    #TODO: Jan!
    #Apply protein length
    pmin = 0.01
    k1 = 30.0*pmin
    k2 = 30.0*pmin
    k3 = 30.0*pmin
    k4 = 30.0*pmin
    k5 = 30.0*pmin


    #Basal activity of inducible promoter
    #TODO: source, RLA, mammalian
    #TODO: Jan, Miha
    alfapc = 0.085
    alfarm = 0.025
    alfapc = alfarm

    #Binding affinity 
    #Literature data. PIP repressor.
    #Tigges.
    KPC = 3.0
    KRM = 3.0
    
    #Cooperativity
    #TODO: source
    npc = self.ghk[1]
    nrm = self.ghk[2]

    #Translation rates (protein/min)
    #TODO: source
    #TODO: protein length adjust
    t57 = 0.02
    t97 = 0.02
    tpc = 0.02
    trm = 0.02
    tnep = 0.02
    tmct = 0.02


    #Degradation rate (mRNA/min)
    #TODO: source
    kdm = 0.0173

    #Degradation rate (protein/min)
    #TODO: source
    kdp = 0.0058

    
    #Degradation rate (inducer/min)
    kdpc = 2/(exp(self.gpc_half_life))
    kdrm = 2/(exp(self.grm_half_life))
    
    #Concentration scaling factor for copies/cell
    #TODO: source
    fv = 1.85e-3*0.5


    #Temporary switch
    #TODO: sources
    r1 = k1*fv
    r2 = k2*fv
    r3 = k3*fv
    r4 = k4*fv
    r5 = k5*fv

  
    '''
    OBSERVED SPECIES
    '''
    HACTm = array([0.0 for i in TIME]) #rapamicin binding protein mRNA
    HACTp = array([0.0 for i in TIME]) #rapamicin binding protein
    HACTa = array([0.0 for i in TIME]) #rapamicin binding protein (active)
    HACTad = array([0.0 for i in TIME]) #rapamicin binding protein (active/draw)

    PBSm = array([0.0 for i in TIME]) #pristinamicin binding protein mRNA
    PBSp = array([0.0 for i in TIME]) #pristinamicin binding protein
    PBSa = array([0.0 for i in TIME]) #pristinamicin binding protein (active)
    PBSad = array([0.0 for i in TIME]) #pristinamicin binding protein (active/draw)

    PC = array([0.0 for i in TIME]) #pristinamicin
    RM = array([0.0 for i in TIME]) #rapamicin

    if self.gpc_half_life > 0:
      PC[self.gpc_start] = 1500
    else:
      l = len(self.gpc_start)
      for i in range(l):
        PC[self.gpc_start[i]:(self.gpc_start[i]+self.gpc_duration[i])] = 1500
        
    if self.grm_half_life > 0:
      RM[self.grm_start] = 1500
    else:
      l = len(self.grm_start)
      for i in range(l):
        RM[self.grm_start[i]:self.grm_start[i]+self.grm_duration[i]] = 1500
        
            
    T97KRABm = array([0.0 for i in TIME])
    T97KRABp = array([0.0 for i in TIME])
    T57KRABm = array([0.0 for i in TIME])
    T57KRABp = array([0.0 for i in TIME])
    
    T97VP16m = array([0.0 for i in TIME])
    T97VP16p = array([0.0 for i in TIME])
    T57VP16m = array([0.0 for i in TIME])
    T57VP16p = array([0.0 for i in TIME])
    
    NEPm = array([0.0 for i in TIME])
    NEPp = array([0.0 for i in TIME])
    MCTm = array([0.0 for i in TIME])
    MCTp = array([0.0 for i in TIME])

    BIFNEPp = []
    BIFMCTp = []
    BIFNEPt = []
    BIFMCTt = []
    
    '''
    SIMULATION
    '''

    
    for t in TIME:
      #Skip first and last step
      if t == 0 or t==self.gTIME-1:
        continue

      dPBSm = G5 * r5 - kdm * PBSm[t-1]
      dPBSp = tpc * PBSm[t-1] - kdp * PBSp[t-1]
      dPBSa = PBSp[t-1] * (PC[t-1]/(Kpc + PC[t-1]))
      dPBSad = 100*PBSp[t-1] * (1 - PC[t-1]/(Kpc + PC[t-1]))

      dHACTm = G5 * r5 - kdm * HACTm[t-1]
      dHACTp = trm * HACTm[t-1] - kdp * HACTp[t-1]
      dHACTa = HACTp[t-1] * (RM[t-1]/(Krm + RM[t-1]))
      dHACTa = HACTp[t-1] * (RM[t-1]/(Krm + RM[t-1]))
      dHACTad = 100*HACTp[t-1] * (RM[t-1]/(Krm + RM[t-1]))
      
      #r3 = fv * k3 * (alfapc + (1-alfapc)*(1-PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc)))
      #r3 = fv * k3 * (alfapc + (alfapc)*(PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc)))
      r3 = fv * k3 * (alfapc + aMIN*(1-alfapc)*(PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc)))   
      r4 = fv * k4 * (alfarm + aMIN*(1-alfarm)*(HACTa[t-1]**nrm/(KRM + HACTa[t-1]**nrm)))
         

      #TAL97:KRAB
      dT97KRABm = G1 * r1 * 1/(1+(T57KRABp[t-1])**hk57/(k57*t57max)**hk57) 
      if T57KRABp[t-1]<self.gthresh:
        dT97KRABm *= a57*(T57VP16p[t-1])**ha57/(Ka57 + T57VP16p[t-1])**ha57
      T97KRABm[t] = T97KRABm[t-1] + dT97KRABm
      T97KRABm[t] = int(T97KRABm[t] > self.gZERO) * T97KRABm[t]
      dT97KRABm = G3*r3 - kdm * T97KRABm[t-1]
      T97KRABm[t] = T97KRABm[t] + dT97KRABm
      dT97KRABp = t97 * T97KRABm[t-1] - kdp * T97KRABp[t-1]

      
      #TAL57:KRAB
      dT57KRABm = G2 * r2 * 1/(1+(T97KRABp[t-1])**hk97/(k97*t97max)**hk97)
      if T97KRABp[t-1]<self.gthresh:
        dT57KRABm *= a97 * (T97VP16p[t-1])**ha97/(Ka97 + T97VP16p[t-1])**ha97
      T57KRABm[t] = T57KRABm[t-1] + dT57KRABm
      T57KRABm[t] = int(T57KRABm[t] > self.gZERO) * T57KRABm[t]
      dT57KRABm = G4*r4 - kdm * T57KRABm[t-1]
      T57KRABm[t] = T57KRABm[t] + dT57KRABm
      dT57KRABp = t57 * T57KRABm[t-1] - kdp * T57KRABp[t-1]
      
      
      #TAL57:VP16
      dT57VP16m   = G1 * r1 * 1/(1+(T57KRABp[t-1])**hk57/(k57*t57max)**hk57)
      if T57KRABp[t-1]<self.gthresh:
        dT57VP16m *= a57*(T57VP16p[t-1])**ha57/(Ka57 + T57VP16p[t-1])**ha57
      T57VP16m[t] = T57VP16m[t-1] + dT57VP16m
      T57VP16m[t] = int(T57VP16m[t] > self.gZERO) * T57VP16m[t]
      dT57VP16m   = G3*r3 - kdm * T57VP16m[t-1]
      T57VP16m[t] = T57VP16m[t] + dT57VP16m
      dT57VP16p   = t57 * T57VP16m[t-1] - kdp * T57VP16p[t-1]

      
      #TAL97:VP16
      dT97VP16m   = G2 * r2 * 1/(1+(T97KRABp[t-1])**hk97/(k97*t97max)**hk97)
      if T97KRABp[t-1]<self.gthresh:
        dT97VP16m  *= a97 * (T97VP16p[t-1])**ha97/(Ka97 + T97VP16p[t-1])**ha97
      T97VP16m[t] = T97VP16m[t-1] + dT97VP16m
      T97VP16m[t] = int(T97VP16m[t] > self.gZERO) * T97VP16m[t]
      dT97VP16m   = G4*r4 - kdm * T97VP16m[t-1]
      T97VP16m[t] = T97VP16m[t] + dT97VP16m
      dT97VP16p   = t97 * T97VP16m[t-1] - kdp * T97VP16p[t-1]
      
      
      #NEPTUN
      dNEPm = G1 * r1 * 1/(1+(T57KRABp[t-1])**hk57/(k57*t57max)**hk57)
      if T57KRABp[t-1]<self.gthresh:
        dNEPm *= a57*(T57VP16p[t-1])**ha57/(Ka57 + T57VP16p[t-1])**ha57
        pass
      dNEPm = dNEPm - kdm * NEPm[t-1]
      dNEPp = tnep * NEPm[t-1] - kdp * NEPp[t-1]

      
      #MCITRIN
      dMCTm = G2 * r2 * 1/(1+(T97KRABp[t-1])**hk97/(k97*t97max)**hk97) 
      if T97KRABp[t-1]<self.gthresh:
        dMCTm  *= a97 * (T97VP16p[t-1])**ha97/(Ka97 + T97VP16p[t-1])**ha97
        pass
      dMCTm = dMCTm - kdm * MCTm[t-1]
      dMCTp = tmct * MCTm[t-1] - kdp * MCTp[t-1]
      

      #PRISTINAMICYN
      dPC = 0 
      if self.gpc_half_life > 0: 
        dPC = -PC[t-1] * kdpc
        
      #RAPAMICYN
      dRM = 0
      if self.grm_half_life > 0: 
        dRM = -RM[t-1] * kdrm
      
      #BIFURCATION
      #All possible phases with dNEPp/dt = 0, dMCTp/dt = 0
      if math.fabs(dNEPp) < self.gZERO and math.fabs(dMCTp) < self.gZERO :
        #Avoid initial state
        if MCTp[t-1] > self.gZERO*1e2 or NEPp[t-1] > self.gZERO*1e2:
          BIFNEPp.append(NEPp[t-1])
          BIFMCTp.append(MCTp[t-1])
      
      
      PBSm[t] = PBSm[t-1] + dPBSm
      PBSp[t] = PBSp[t-1] + dPBSp
      HACTm[t] = HACTm[t-1] + dHACTm
      HACTp[t] = HACTp[t-1] + dHACTp   
      PBSa[t] = PBSa[t] + dPBSa
      PBSad[t] = PBSad[t] + dPBSad
      HACTa[t] = HACTa[t] + dHACTa
      HACTad[t] = HACTad[t] + dHACTad
      T97KRABp[t] = T97KRABp[t-1] + dT97KRABp
      T57KRABp[t] = T57KRABp[t-1] + dT57KRABp
      T97VP16p[t] = T97VP16p[t-1] + dT97VP16p
      T57VP16p[t] = T57VP16p[t-1] + dT57VP16p
      NEPm[t] = NEPm[t-1] + dNEPm
      NEPp[t] = NEPp[t-1] + dNEPp
      MCTm[t] = MCTm[t-1] + dMCTm
      MCTp[t] = MCTp[t-1] + dMCTp  
      
      if self.gpc_half_life > 0 and not t in self.gpc_start:
        PC[t] = PC[t-1] + dPC
      if self.grm_half_life > 0 and not t in self.grm_start:
        RM[t] = RM[t-1] + dRM

      #Preserve above zero
      PBSm[t] = int(PBSm[t] > self.gZERO) * PBSm[t] + int(PBSm[t] < self.gZERO) * self.gZERO
      PBSp[t] = int(PBSp[t] > self.gZERO) * PBSp[t] + int(PBSp[t] < self.gZERO) * self.gZERO
      HACTm[t] = int(HACTm[t] > self.gZERO) * HACTm[t] + int(HACTm[t] < self.gZERO) * self.gZERO
      HACTp[t] = int(HACTp[t] > self.gZERO) * HACTp[t] + int(HACTp[t] < self.gZERO) * self.gZERO
      PBSa[t] = int(PBSa[t] > self.gZERO) * PBSa[t] + int(PBSa[t] < self.gZERO) * self.gZERO
      PBSad[t] = int(PBSad[t] > self.gZERO) * PBSad[t] + int(PBSad[t] < self.gZERO) * self.gZERO
      HACTa[t] = int(HACTa[t] > self.gZERO) * HACTa[t] + int(HACTa[t] < self.gZERO) * self.gZERO
      HACTad[t] = int(HACTad[t] > self.gZERO) * HACTad[t] + int(HACTad[t] < self.gZERO) * self.gZERO
      T97KRABp[t] = int(T97KRABp[t] > self.gZERO) * T97KRABp[t] + int(T97KRABp[t] < self.gZERO) * self.gZERO
      T57KRABp[t] = int(T57KRABp[t] > self.gZERO) * T57KRABp[t] + int(T57KRABp[t] < self.gZERO) * self.gZERO
      T97VP16p[t] = int(T97VP16p[t] > self.gZERO) * T97VP16p[t] + int(T97VP16p[t] < self.gZERO) * self.gZERO
      T57VP16p[t] = int(T57VP16p[t] > self.gZERO) * T57VP16p[t] + int(T57VP16p[t] < self.gZERO) * self.gZERO
      NEPm[t] = int(NEPm[t] > self.gZERO) * NEPm[t] + int(NEPm[t] < self.gZERO) * self.gZERO
      NEPp[t] = int(NEPp[t] > self.gZERO) * NEPp[t] + int(NEPp[t] < self.gZERO) * self.gZERO
      MCTm[t] = int(MCTm[t] > self.gZERO) * MCTm[t] + int(MCTm[t] < self.gZERO) * self.gZERO
      MCTp[t] = int(MCTp[t] > self.gZERO) * MCTp[t] + int(MCTp[t] < self.gZERO) * self.gZERO
      PC[t] = int(PC[t] > self.gZERO) * PC[t] + int(PC[t] < self.gZERO) * self.gZERO
      RM[t] = int(RM[t] > self.gZERO) * RM[t] + int(RM[t] < self.gZERO) * self.gZERO

    if self.gDEBUG:
      #return self.bifurcationAnalysis(BIFNEPp,BIFMCTp, NEPp, MCTp)
      pass
    return self.plotResults(TIME, T97KRABm, T97KRABp, T57KRABm, T57KRABp, T97VP16m, T97VP16p, T57VP16m, T57VP16p,  PBSm, PBSp, PBSad, HACTm, HACTp, HACTad, PC, RM, NEPm, NEPp, MCTm, MCTp, BIFNEPp, BIFNEPt, BIFMCTp, BIFMCTt)
    


      

